#ifndef STEIN_GCD_H
#define STEIN_GCD_H

template<typename T>
T stein_gcd_recursive(T a, T b)
{
	// simple cases (termination)
	if (a == b)
		return a;
	if (a == 0)
		return a;
	if (b == 0)
		return b;

	// look for factors of 2
	if (~a & 0x01) { // a is even
		if (b & 0x01) // b is odd
			return stein_gcd_recursive(a >> 1, b);
		else // both a and b are even
			return stein_gcd_recursive(a >> 1, b >> 1) << 1;
	}
	else if (~b & 0x01) {// a is odd, b is even
		return stein_gcd_recursive(a, b >> 1);
	}
	else { //a is odd, b is odd
		if (a > b)
			return stein_gcd_recursive((a - b) >> 1, b);
		else
			return stein_gcd_recursive((b - a) >> 1, a);
	}
}

#endif